##########################################################################################

library(dplyr)
library(ggplot2)
library(data.table)
library(RColorBrewer)
library(optparse)
library(ggpubr)
library(ggsci)

##########################################################################################

option_list <- list(
    make_option(c("--input_file"), type = "character") ,
    make_option(c("--images_path"), type = "character")
)

if(1!=1){
    
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"
    input_file <- paste(work_dir,"/images/mutBurden/mutBurden.tsv",sep="")
	images_path <- paste(work_dir,"/images/mutBurden",sep="")

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

input_file <- opt$input_file
images_path <- opt$images_path

dir.create(images_path , recursive = T)

###########################################################################################

col <- c(
  brewer.pal(9,"YlGnBu")[6],
  rgb(234,106,79,alpha=255,maxColorValue=255),
  rgb(203,24,30,alpha=255,maxColorValue=255),
  rgb(255,0,0,alpha=255,maxColorValue=255)
  )

names(col) <- c("IM" , "IGC" , "DGC" , "GC")
col <- col[1:3]
class_order <- c("IGC" , "DGC" , "IM")

Variant_Types <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")

###########################################################################################

dat_input <- data.frame(fread( input_file ))

###########################################################################################
## 多个样本负荷取中位数
dat_plot2 <- dat_input %>%
group_by( Patient , Age , Gender , Tobacco , Alcohol , Alcohol_frequncy , PickleFood , HP , Class , Type , MS_Type ) %>%
summarize( BurdenAll = median(BurdenAll) ,BurdenExon = median(BurdenExon) )

dat_plot2$Class <- factor( dat_plot2$Class , levels = class_order , order = T )

###########################################################################################
plotBurden <- function(dat_plot_tmp = dat_plot_tmp , out_name = out_name , my_comparisons_1 = my_comparisons_1 , width= width){

    plot <- ggplot( dat_plot_tmp , aes( x = useCol , y = MutBurden_use , color = useCol ) ) +
        geom_boxplot(alpha =1 , outlier.size=0 , size = 1.5 , width = 0.6) +
        geom_jitter(position = position_jitter(0.17) , size = 1 , alpha = 0.7) +
        facet_grid(.~Class,space='free_x',scales='free_x') +
        scale_color_npg() +
        xlab(NULL) +
        ylab("Mutation rate per MB")+
        theme_bw() +
        stat_compare_means(comparisons = my_comparisons_1,method = "wilcox.test") +
        theme(panel.background = element_blank(),#设置背影为白色#清除网格线
            legend.position ='none',
            legend.title = element_blank() ,
            panel.grid.major=element_line(colour=NA),
            plot.title = element_text(size = 8,color="black",face='bold'),
            legend.text = element_text(size = 10,color="black",face='bold'),
            axis.text.y = element_text(size = 10,color="black",face='bold'),
            axis.title.x = element_text(size = 10,color="black",face='bold'),
            axis.title.y = element_text(size = 10,color="black",face='bold'),
            axis.ticks.x = element_blank(),
            axis.text.x = element_text(size = 10,color="black",face='bold',angle = 45, hjust = 1) ,
            axis.line = element_line(size = 0.5)) 

     
    ggsave(file=out_name,plot=plot,width=width,height=5)

}

baseuse <- "Alcohol_frequncy"
mutType <- "cds"

dat_plot_tmp <- data.frame(dat_plot2)
dat_plot_tmp$useCol <- dat_plot2[[baseuse]]
dat_plot_tmp <- subset( dat_plot_tmp , Class == "IM" | Type != "IM + IGC + DGC" )
dat_plot_tmp <- subset( dat_plot_tmp , MS_Type != "MSI" | Class == "IM" )
dat_plot_tmp <- subset( dat_plot_tmp , !is.na(useCol) )
dat_plot_tmp$useCol <- ifelse( dat_plot_tmp$useCol == "A" , "Never" , dat_plot_tmp$useCol )
dat_plot_tmp$useCol <- ifelse( dat_plot_tmp$useCol == "B" , "Sometime" , dat_plot_tmp$useCol )
dat_plot_tmp$useCol <- ifelse( dat_plot_tmp$useCol == "C" , "<4/month" , dat_plot_tmp$useCol )
dat_plot_tmp$useCol <- ifelse( dat_plot_tmp$useCol == "D" , "1/week" , dat_plot_tmp$useCol )
dat_plot_tmp$useCol <- ifelse( dat_plot_tmp$useCol == "E" , "2-3/week" , dat_plot_tmp$useCol )
dat_plot_tmp$useCol <- ifelse( dat_plot_tmp$useCol == "F" , "4-5/week" , dat_plot_tmp$useCol )
dat_plot_tmp$useCol <- ifelse( dat_plot_tmp$useCol == "G" , "EveryDay" , dat_plot_tmp$useCol )
dat_plot_tmp$useCol <- factor( dat_plot_tmp$useCol , 
    levels = c("Never" , "Sometime" , "<4/month" , "1/week" , "2-3/week" , "4-5/week" , "EveryDay") , order = T )

my_comparisons_1 <- list( 
    c(1, 2) , c(1,4) , c(1,5) ,
    c(2,4) , c(2,5) ,
    c(4,5) , c(6,5)  
    )

dat_plot_tmp$MutBurden_use <- dat_plot_tmp$BurdenExon

table(dat_plot_tmp$useCol)

for(class in c("IM" , "IGC" , "DGC")){

    if(class=="DGC"){
        my_comparisons_1 <- list( 
        c(1, 2) , c(1,4) , c(1,5) ,
        c(2,4) , c(2,5) ,
        c(4,5) 
        )
    }
    dat_plot_tmp_use <- subset( dat_plot_tmp , Class == class )
    out_name <- paste0( images_path , "/mutBurden.",baseuse,".",mutType,".",class,".pdf" ) 
    plotBurden(dat_plot_tmp = dat_plot_tmp_use , out_name = out_name , my_comparisons_1 = my_comparisons_1 , width=4)
}

###########################################################################################
## IGC的IM和DGC的IM是否存在差异
baseuse <- "Alcohol_frequncy"
mutType <- "cds"

dat_plot_tmp <- data.frame(dat_plot2)
dat_plot_tmp$useCol <- dat_plot2[[baseuse]]
dat_plot_tmp <- subset( dat_plot_tmp , !is.na(useCol) )
dat_plot_tmp1 <- dat_plot_tmp

## 比较IGC的IM和DGC的IM
dat_plot_tmp <- subset( dat_plot_tmp1 , Class == "IM" & Type != "IM + IGC + DGC" )
dat_plot_tmp$Class <- as.character(dat_plot_tmp$Class)
dat_plot_tmp[dat_plot_tmp$Type=="IM + IGC" , "Class"] <- "IM(IGC)"
dat_plot_tmp[dat_plot_tmp$Type=="IM + DGC" , "Class"] <- "IM(DGC)"

dat_plot_tmp2 <- subset( dat_plot_tmp1 , Class == "IM" & Type == "IM + IGC + DGC" )
dat_plot_tmp2$Class <- "IM(IGC_DGC)"
dat_plot_tmp <- rbind( dat_plot_tmp2 , dat_plot_tmp )

if(baseuse=="HP"){
    dat_plot_tmp[dat_plot_tmp$useCol=="Negative","useCol"] <- "HP_N"
    dat_plot_tmp[dat_plot_tmp$useCol=="Positive","useCol"] <- "HP_Y"
    dat_plot_tmp$useCol <- factor( dat_plot_tmp$useCol , levels = c("HP_Y" , "HP_N") , order = T )
}else if(baseuse=="Tobacco"){
    dat_plot_tmp$useCol <- factor( dat_plot_tmp$useCol , levels = c("Smoke" , "No") , order = T )
}

dat_plot_tmp$useCol <- ifelse( dat_plot_tmp$useCol == "A" , "Never" , dat_plot_tmp$useCol )
dat_plot_tmp$useCol <- ifelse( dat_plot_tmp$useCol == "B" , "Sometime" , dat_plot_tmp$useCol )
dat_plot_tmp$useCol <- ifelse( dat_plot_tmp$useCol == "C" , "<4/month" , dat_plot_tmp$useCol )
dat_plot_tmp$useCol <- ifelse( dat_plot_tmp$useCol == "D" , "1/week" , dat_plot_tmp$useCol )
dat_plot_tmp$useCol <- ifelse( dat_plot_tmp$useCol == "E" , "2-3/week" , dat_plot_tmp$useCol )
dat_plot_tmp$useCol <- ifelse( dat_plot_tmp$useCol == "F" , "4-5/week" , dat_plot_tmp$useCol )
dat_plot_tmp$useCol <- ifelse( dat_plot_tmp$useCol == "G" , "EveryDay" , dat_plot_tmp$useCol )
dat_plot_tmp$useCol <- factor( dat_plot_tmp$useCol , 
    levels = c("Never" , "Sometime" , "<4/month" , "1/week" , "2-3/week" , "4-5/week" , "EveryDay") , order = T )

dat_plot_tmp$Class <- factor( dat_plot_tmp$Class , levels = c("IM(IGC)" , "IM(DGC)","IM(IGC_DGC)") , order = T )
dat_plot_tmp$MutBurden_use <- dat_plot_tmp$BurdenExon

my_comparisons_1 <- list( 
    c(1, 2) , c(1,4) , c(1,5) ,
    c(2,4) , c(2,5) ,
    c(4,5) , c(6,5)  
    )

for(class in c("IM(IGC)" , "IM(DGC)")){
    if(class=="IM(DGC)"){
        my_comparisons_1 <- list( 
        c(1, 2) , c(1,4) , c(1,5) ,
        c(2,4) , c(2,5) ,
        c(4,5) 
        )
    }
    dat_plot_tmp_use <- subset( dat_plot_tmp , Class == class )
    out_name <- paste0( images_path , "/mutBurden.",baseuse,".",mutType,".",class,".divide.pdf" ) 
    plotBurden(dat_plot_tmp = dat_plot_tmp_use , out_name = out_name , my_comparisons_1 = my_comparisons_1 , width=4)
}
